/**
 * @file Polynomial.H
 * @brief 一个描述一元实系数多项式的模板类，可用于加，减，乘，求导，求值，与格式化输出。
 * T 可以是 数， 也可以是 方阵 等具有域与线性空间结构的抽象类。
 * @author htli (3180102114@zju.edu.cn)
 * @version 1.0
 * @date 2021-10-22
 * 
 * @copyright Copyright (c) 2021  linhuo2020
 * 
 */
#ifndef _PROJECT2_POLYNOMIAL_H_
#define _PROJECT2_POLYNOMIAL_H_

#include <iostream>
#include <vector>
#include <cmath>

#define Template template <class T>
#define Poly Polynomial<T>

#ifndef _TYPEDEF_REAL_
#define _TYPEDEF_REAL_
///根据需要更改 double 为 long double 或 float
typedef double Real;
typedef std::vector<Real> RealVect;
typedef std::vector<RealVect> RealBox;
#else
//do nothing
#endif

Template class Polynomial;
Template const Poly operator+(const Poly &_p1, const Poly &_p2);
Template const Poly operator-(const Poly &_p1, const Poly &_p2);
Template const Poly operator*(const Poly &_p1, const Poly &_p2);
Template const Poly operator*(const Real &_x, const Poly &_p1);
Template std::ostream &operator<<(std::ostream &os, const Poly &_p);

Template class Polynomial
{
protected:
    ///依次存储常数项至最高次项的实系数（末尾项若不是常数项，系数必须非0）
    ///默认创建零函数
    RealVect coef_ = {0};
    ///系数的精度（两数之差小于容忍度时会被认为一个数）
    Real tolerance_ = 1e-7;

public:
    Polynomial(){};
    ///创建时给定系数
    Polynomial(const RealVect _coef);
    ~Polynomial(){};
    /// 重改系数
    void ResetCoef(const RealVect _coef);
    /// 简化系数以确保最高项非0
    void SimplifyCoef();
    ///获取阶数
    const unsigned int getDeg() const { return (coef_.size() - 1); }
    ///获取系数
    const RealVect getCoef() const { return coef_; };
    ///获取系数精度
    const Real getTol() const { return tolerance_; };
    ///计算多项式作用于_x上的结果
    const T eval(const T &_x) const;
    ///计算多项式相加的结果
    friend const Poly operator+<>(const Poly &_p1, const Poly &_p2);
    ///计算多项式相减的结果
    friend const Poly operator-<>(const Poly &_p1, const Poly &_p2);
    ///计算多项式相乘的结果
    friend const Poly operator*<>(const Poly &_p1, const Poly &_p2);
    ///计算多项式数乘的结果
    friend const Poly operator*<>(const Real &_x, const Poly &_p1);
    ///格式化输出多项式
    friend std::ostream &operator<<<>(std::ostream &os, const Poly &_p);
    ///计算多项式的一阶导数
    const Poly derivate() const;
};

Template Poly::Polynomial(const RealVect _coef)
{
    this->ResetCoef(_coef);
};

Template void Poly::ResetCoef(const RealVect _coef)
{
    coef_ = _coef;
    this->SimplifyCoef();
};

Template void Poly::SimplifyCoef()
{
    while (coef_.size() > 1 && fabs(coef_.back()) < tolerance_)
        coef_.pop_back();
};

Template const T Poly::eval(const T &_x) const
{
    //寄。这下没法对方阵用多项式了。(或许这么定义还有抢救的余地？)
    T pow_x = (T)(1);
    T px = (T)(0);
    const int deg = this->getDeg();
    for (int i = 0; i <= deg; i++)
    {
        px += coef_[i] * pow_x;
        pow_x *= _x;
    };
    return px;
};

Template const Poly Poly::derivate() const
{
    RealVect dPCoef = this->getCoef();
    const int dPDeg = this->getDeg();
    for (int i = 0; i < dPDeg; i++)
        dPCoef[i] = (i + 1) * dPCoef[i + 1];
    dPCoef[dPDeg] = 0;
    Poly p(dPCoef);
    return p;
};

Template const Poly operator+(const Poly &_p1, const Poly &_p2)
{
    RealVect p1Coef = _p1.getCoef();
    RealVect p2Coef = _p2.getCoef();
    const int M = std::max(_p1.getDeg(), _p2.getDeg()) + 1;
    p1Coef.resize(M, 0);
    p2Coef.resize(M, 0);
    RealVect pCoef(M);
    for (int i = 0; i < M; i++)
        pCoef[i] = p1Coef[i] + p2Coef[i];
    Poly p(pCoef);
    return p;
};

Template const Poly operator-(const Poly &_p1, const Poly &_p2)
{
    Poly p3 = -1 * _p2;
    Poly p = _p1 + p3;
    return p;
};

Template const Poly operator*(const Poly &_p1, const Poly &_p2)
{
    RealVect p1Coef = _p1.getCoef();
    RealVect p2Coef = _p2.getCoef();
    RealVect pCoef;
    const int M = (_p1.getDeg()) + (_p2.getDeg()) + 1;
    p1Coef.resize(M, 0);
    p2Coef.resize(M, 0);
    pCoef.resize(M, 0);
    for (int i = 0; i < M; i++)
        for (int j = 0; j < M - i; j++)
            pCoef[i + j] += p1Coef[i] * p2Coef[j];
    Poly p(pCoef);
    return p;
};

Template const Poly operator*(const Real &_x, const Poly &_p1)
{
    RealVect p2Coef(1, _x);
    Poly p2(p2Coef);
    Poly p = _p1 * p2;
    return p;
};

Template std::ostream &operator<<(std::ostream &os, const Poly &_p)
{
    const RealVect pCoef = _p.getCoef();
    const int pDeg = pCoef.size() - 1;
    const Real tol = _p.getTol();
    os << "px = ";
    os << pCoef[pDeg] << "*x.^" << pDeg;
    for (int i = pDeg - 1; i >= 0; i--)
    {
        if (pCoef[i] > tol)
        {
            os << " +" << pCoef[i] << "*x.^" << i;
        }
        else if (pCoef[i] < (-tol))
        {
            os << " " << pCoef[i] << "*x.^" << i;
        };
    }
    os << " ";
    return os;
}

#else
//do nothing
#endif
